The Kasteleyn transition in three dimensions: spin ice in a [100] field 
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We discuss the nearest neighbour spin ice model in the presence of a magnetic field placed along 
the cubic [100] direction. As recently shown in Phys. Rev. Lett. 100, 067207, 2008, the symmetry 
sustaining ordering transition observed at low temperature is a three dimensional Kasteleyn transi- 
tion. We confirm this with numerical data using a non-local algorithm that conserves the topological 
constraints at low temperature and from analytic calculations from a Bethe lattice of corner sharing 
tetrahedra . We present a thermodynamic description of the Kasteleyn transition and discuss the 
relevance of our results to recent neutron scattering experiments on spin ice materials. 
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This transition was introduced by Kasteleyn [l[ in the context of ordering of hard core 2D dimcrs lying on the bonds 
\q [ of a honeycomb lattice. Applying local chemical potentials, /ij, for dimers oriented in one of the i = 1,2,3 possible 
directions leads to a singular ordering transition for finite values of the fugacities Zi = exp{— j3 Hi) . The transition 
can also occur for dimers on three dimensional equivalents of the honeycomb lattice such as the brick lattice [2], 
while the exponents can be determined in dimension d through mapping onto a directed polymer problem Q. It has 
been extensively used to describe trans-gauche structural transitions in polymerized lipid bilayers and the resulting 
theory provides a good qualitative description of the singular density change occurring in solvent bilayer systems as 
I 1 a function of temperature A Kasteleyn transition also occurs in the five vertex model on a square lattice [f|, 
which maps onto the square ice model, satisfying the ice rules in two dimensions with direct relevance to spin ice 
materials [6J. It has been shown that the same Kasteleyn transition exists in the nearest neighbour model for spin 
ice materials on a pyrochlore lattice, with magnetic field directed close to the body centred cubic [111] direction [3] ■ 
Placing the field along this direction isolates layers orthogonal to the field, each of which forms a kagome lattice 
of directed tetrahedra whose vertices lie alternately above and below the plane. The magnetic field fixes only one 
I , out of the four sublattices with the result that the extensive ground state degeneracy of zero field is only partially 
lifted. The ground state entropy remains finite and the kagome spin problem maps exactly onto Kasteleyn's original 
dimer problem on the Honeycomb lattice. Tilting the field off the [111] direction, giving a finite perpendicular field 
component H±, is equivalent to applying chemical potentials /ij to the dimers and the spins order in the kagome 
planes for a finite value of H^/ksT via Kasteleyn's transition. Ground state configurations of this model require 
the imposition of local constraints, which lead to power law dipole like correlations [&, 9] . Recent neutron studies of 
Ho2Ti27 with this field configuration (To| have shown experimentally the existence of these correlations in the form of 
l^-j ■ well defined pinch points in reciprocal space for the intensity of the elastic diffuse scattering. In two dimensions these 
ly" - ) \ points should include logarithmically divergent intensities. In the presence of a perpendicular field component the 
t-H ■ singularities are predicted 0] to drift with changing field strength, H±. It has recently been shown that the ordering 
transition seen when a field is placed along the [100] direction fll| is a three dimensional Kasteleyn transition [l2l ]. 
£T} • Magnetic measurements on Holmium and Dysprosium Titanate[13Lll4| indicate that experimental systems could show 
strong features of the transition, but as we will see below, accessing it while remaining in thermodynamic equilibrium 
is a considerable experimental challenge. 
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In this paper we review the behaviour of the nearest neighbour spin ice model, with magnetic field placed along the 
cubic [100] direction. Nearest neighbour moments (spins) of unit length are constrained to lie along the body centred 
crystal field directions of the tetrahedra: Si = ±di and interact with a ferromagnetic exchange interaction J and with 
an external magnetic field, H in reduced units. The Hamiltonian Q, is 

{ij) i 



II. THE KASTELEYN TRANSITION IN SPIN ICE IN A [100] FIELD 



The crystal fields frustrate the ferromagnetic exchange leading to a highly degenerate ground state satisfying the so 
called ice rules [la] with two spins pointing in and two spins pointing out of each tetrahedron. This is the topological 
constraint which takes the form V.S'(r) = Q in a continuum description. In zero field there is no phase transition 
down to zero temperature but, in the presence of a magnetic field and in the limit where the topological constraints 
are rigorously applied the system orders at low temperature via a three dimensional [l2j Kasteleyn transition. 

The remarkable thing about the Kasteleyn transition is that it occurs when the system is rigorously constrained to 
the manifold of spin ice states such that the internal energy for each microstate in the reduced phase space has the 
same value, U — —NJ/3. This allows for a particularly simple thermodynamic description which we present below. 
Just as for a paramagnet the magnetic Helmholtz free energy for the constrained system has an entropic component 
only: F(N,T, M) = —TS(N,T,M), where M is the magnetic order parameter. In a regular paramagnetic system S 
goes to zero, as M reaches its maximum value, M max but it does so with infinite slope and hence an infinite value for 
H/ksT, where H = \H\ is the field conjugate to M: 

H- 1 dF (2) 
H -~NdM- (2) 

That is, there is no phase transition. However the hidden divergence free constraint can drive the entropy to zero for 
finite H/ksT giving a singular longitudinal susceptibility \ — dM/dH-this is the Kasteleyn transition. The Gibbs 
potential is defined G* = F — NMH from which one can construct a Landau free energy, Gl by expanding G* in 
small x = M max — M near the transition: 

^ = (H- H K )x + ^-Tx 2 + ^T.t 3 + • • • - M max h, (3) 

where Hx is the critical field for T fixed and where ai are constants. If «2 > one finds the mean field critical 
behaviour for fixed H: 

x = x (T-T K f, (3 = 1, (4) 

where Tk is the transition temperature for fixed H. One also finds that the specific heat, Ch, and susceptibility, 
X, jump discontinuously to zero at Tk. For T < Tk G*(M) has no minimum for < M < M max . dG*/dM is 
finite and negative at M = M maxi consistently with there being zero fluctuations and the entropy being strictly zero: 
the system is perfectly ordered and represents a complete vacuum for excitations. Hence the transition is violently 
asymmetric, with fluctuations above the transition but not below. Through this asymmetry the transition has often 
been described as having both I s * and 2 nd order characteristics[l6|. However from the above thermodynamics it is 



clearly second order, or continuous, within the Ehrenfest classification[l7|. 

Zero entropy can only be maintained at finite temperature if the topological constraints are rigorously applied 
(J >> ksT) so that the lowest energy excitations off the perfectly ordered state are extended, with infinitely large 
energy in the thermodynamic limit. To see that this is the case consider a strong field at low temperatures which moves 



the system from a disordered spin ice state to the long range ordered state illustrated by the black spins in Fig 1(a) 
Relaxing back to a disordered state while maintaining the ice rules requires the flipping of a string of spins spanning 
the entire system (in the following discussion we impose periodic boundary conditions, in which case the string closes 
through the boundaries). To see this consider flipping a single spin out of the ordered state. This breaks the ice rules 
creating a pair of "topological defects" (la |. tetrahedra with three spins in and one out and vise versa, which can 
only be annealed by interacting with each other. Creating the defects costs magnetic enthalpy 4 J/3 + 2_ff/v / 3, but 
moving them further apart costs no further exchange energy in the nearest neighbour model. The exchange energy 
can be regained if the ensuing string passes right through the system and the topological defects destroy each other, 
returning the system to the constrained manifold of spin ice states. Such a string is illustrated by the red spins in 



Figure 1 (a) Ignoring the exchange energy, each step of the string costs Zeeman energy 2H/ ^/(3) but also represents 



an entropy gain of ln(2) as the string can leave the tetrahedron in one of two directions. If we define L as the number 
of spins in the [100] direction, the total free enthalpy change for placing a string in the system is then 

5G = L (^L -Tln2 

oc L (Tk — T) , where 

_ 2H/y/{3) 

TK In 2 ' (5) 

and strings appear above the temperature Tk only. The transition could be first order if strings attracted each other 
but two strings passing though the same tetrahedron clearly have a reduced entropy over that for two separated 
strings. Hence, strings repel each other and the result is a continuous Kasteleyn transition at temperature Tk- 




FIG. 1: (a) The pyrochlore lattice with spins satisfying the ice rules. H is given along the [100] direction. Black spins: long 
range ordered state, red spins: a string of nipped spins. 

Figure 1 (b) (Color online) Magnetization per spin M vs. T/Tk obtained from simulations for the pyrochlore lattice 
(dots) and analytically on the Bethe lattice (solid lines) for H/J = 10 _3 (»), 0.13(a) and 0.58(H). 



A consequence of the topological constraints and the broken magnetic symmetry is that a mapping exists between 
the strings of returned spins and world lines for bosons living in a d — 1 dimensional space perpendicular to the field 
[l9j . The pyrochlore lattice has six fold cubic symmetry, which is reduced by the magnetic field, giving a unique 
definition for the arrow of time of the world lines along an imaginary time axis. The bosonic nature of the particles 
shows up in the passage of two strings through a single tetrahedron . The entropic contribution of the tetrahedron is 
zero and it is impossible to index strings "1" and "2" either side of the encounter without a fixed convention. However, 
with such a convention an exact mapping exists between spin and string configurations. In the absence of a field, 
there are six equivalent definitions for the arrow of time and a further convention is required to define the strings. 
In reference [121 it was shown that through the definition of strings as world lines, the Kasteleyn transition in 
dimension d maps onto a quantum phase transition in dimension d — 1. The applied field corresponds to the boson 
chemical potential, /u *~ (h — hx) and the transition is from a Bose condensed to a vacuum state at zero temperature. 
The upper critical dimension for this transition is d = 3 [20], which is compatible with previous analytic studies of the 
Kasteleyn transition in three dimensional dimer models[2j. Hence, for the three dimensional Kasteleyn transition we 
expect near mean field behaviour, with logarithmic corrections to scaling, while in two dimensions fluctuations should 
be important. This is indeed the case: the analytic solution of the two dimensional problem^ gives a value for the 
exponent, (3 = 1/2. 

We have tested these ideas using a non-local Monte Carlo algorithm which flips connected lines of spins [H, 
l2l| . The lines can form closed loops, either within the sample or through the boundaries, or they can transport 
topological defects non-locally. In the case of closed loops only, the algorithm restricts the system to fully constrained 
spin ice states, allowing an ergodic exploration of this manifold without passing through energetically unfavorable 



unconstrained states. In Figure 1(b) we show M = -k\ £\ Si\ vsT for fixed H. The points represent data generated by 
simulation using the non-local loop algorithm, while solid lines show analytic data from a calculation on a tree structure 
(Bethe lattice) of connected tetrahedron . The lack of closed loops here, allows for an essentially analytic solution 
of the problem [221 ]. For H/J = 0.001 the data show the characteristics of the Kasteleyn transition: asymmetric 
singular behaviour, with the magnetisation arriving at its saturation value M max = 1/V3 precisely at T = 7k- This 



is confirmed in figure 2(a) where the magnetisation difference, x = M max — M and the specific heat (and hence the 
susceptibility) are both seen to be singular. As predicted through the mapping to the quantum phase transition, 
the simulation data for Ch fit accurately to a logarithmic divergence above the transition, while x approaches zero 
linearly with T — Tj< and has logarithmic corrections for larger x: 

x = AM oc t (1 - aln(t)) ^ C h oc - ln(i) (6) 

Detailed calculations for the Bethe lattice confirm the mean field picture developed phenomenologically above [22| . 

Previous numerical work found indications of a 1 st order phase transition at finite field [ll|. The observed first order 
nature is a consequence of the loss of ergodicity due to the single spin flip dynamics used at that time. Interestingly 
experimental results [HI] reflected the same out-of-equilibrium features, as natural dynamics is also local. However, 
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FIG. 2: (a) (Color online) Specific heat Ch and relative magnetisation AM as a function of t obtained from simulations (• for 
Ch and ■ for AM) and analytically (solid line for Ch) for H — 10~ 3 J. The dashed lines are a fit to logarithmic singularities 
given in eq. ©. 

Figure 2 (b) (Color online) Finite-size scaling of spin correlations for H/H c = 0.33(H), 0.66(a) and 0.93(»), on log- 
log scale. Lines have slopes 3 and 1. 



with the introduction of the non-local algorithm the numerical results very clearly confirm that the transition is in 
fact a Kasteleyn transition similar to that observed for dimers on the brick lattice [2|. 

For larger fields the transition is rounded, this is because the topological constraints are no longer rigorously imposed 
and the system can support a finite concentration of topological defects. Above a characteristic scale fixed by their 
separation the system will become regularly paramagnetic and this is just what is observed in the rounding of the 
transition. 

Despite the fact that the constraints induce long range correlations the ferromagnetic susceptibility is not divergent 
away from the Kasteleyn transition. The coupling between spin and configurational space for dipolar interactions 
leads to both direction and sublattice dependence for the amplitude of the spin spin correlation function. This gives 
a finite susceptibility at zero wave vector, = 0, but very strong variation at finite Q, leading to the so called pinch 
points on the Brillouin zone boundaries |23[. Although logarithmically divergent in 2D the intensity of the pinch 
points in three dimensions remains finite in zero field. Applying the appropriate magnetic field in the 2D problem 
breaks the square symmetry of zero field, introducing a length scale related to the symmetry breaking. As a result, the 
pinch points move in towards the Brillouin zone centre. Their arrival at the zone centre corresponds to the Kasteleyn 
transition. In 3D the pinch points remain fixed in reciprocal space but their intensity evolves with applied field, with 
the point at the Brillouin zone centre developing a logarithmic divergence at the transition. This evolution can be 
seen directly in real space through spin-spin correlation functions C{r) —< Si(r).Si(0) >. In Figure [2(b)] we show 
finite size scaling plots of C(z = L/2) for three dimensional systems of size L for magnetic field along the z axis, from 
which we can extract the power law behaviour of C(r) in the thermodynamic limit. For small field we find a 1/r 3 
behaviour consistently with the dipolar correlations leading to the pinch points. As one approaches the transition 
we observe a crossover to 1/r dependence. This is the regime where spin fluctuations along z are due to multiple 
encounters with a single string. Here, mapping either to the quantum problem or to a random walk problem lead to 
the 1/r behaviour observed in the figure [12j, [l9|. 

The loop algorithm allows particularly accurate studies of finite size effects close to the transition where the density 
of strings is small. In the ordered regime, the free energy cost for the excitation of a string @, 5G is dominated by 
the energy E str ing — 2 H L/ V3 hence the Boltzmann weight for string excitation decreases exponentially with system 
size: wi oop cx exp (—(3E str ing) [HJ. As the magnetisation difference x is proportional to the number of strings one 
can observe this exponential finite size scaling by making an accurate measure of x in the transition region. In Fig. 
(3(a) ) we show a; on a logarithmic scale, as a function of Tk/T, for different system sizes. The data, averaged over up 
to 10' configurations fit the exponential scaling remarkably well down to x ~ 10 . This analysis shows convincingly 
that the mechanism for the transition is indeed the generation of extensive excitations developed in this section. 
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FIG. 3: (Color online) Relative magnetisation AM as a function of H/Hk obtained numerically (•) for different sizes of the 
system: L = 20, 40, 100, 200. The straight lines are a fit to the Boltzmann factor e (L ln2) h/Hk . 

Figure 3 (b) M vs B obtained analytically (solid lines) and experimentally (dots) by neutron scattering on the com- 
pound Ho2Ti207 at T — 1.2K. For fitting details, see main text. 



III. DISCUSSION 



As the experimental systems fall out of equilibrium at low temperature, it is difficult to compare theory with 
experiment very close to the Kasteleyn transition. But at reasonably low temperatures, where the concentration of 
topological defects is low but sufficiently high to maintain ergodicity, it is possible to observe a rounded transition 
and to compare data with our theory. As an example we show, in Fig. |3(b)1 one of the [100] field scans performed 
at fixed low temperature in Holmium Titanate [la ]. The data is for T = 1.2K. This temperature is about 2/3 of the 
effective nearest neighbour exchange constant J e ffS 2 /ks ~ 1.8K estimated for Ho2Ti207 [Hj]- The magnetisation 
flattens off abruptly at the saturation value, M exp — IO/ib/V^ in a way reminiscent of the Kasteleyn transition 
but is indeed rounded near saturation. Also shown in the figure is our data calculated analytically from the Bcthe 
lattice approximation for the pyrochlore system with the same ratio of J /ksT '. The scale of the magnetic field (that 
is changing from \H\ in reduced units to B in Tesla) was estimated after demagnetisation effects were taken into 
account on the experimental data. A final scale factor of 1.6 was required to make the comparison. The comparison 
is qualitatively good, suggesting that the data of Ref. l3j do provide evidence for a rounded Kasteleyn transition in 
Ho2Ti207, which could be more closely approached in future experiments. 
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